require(Rmisc) # Package for doing data summaries. We'll use its summarySE function later to get means and SEs of various datasets for plotting.
## Loading required package: Rmisc
## Loading required package: lattice
## Loading required package: plyr
require(sciplot) # The bargraph.CI function in this package makes barplots w/ error bars a breeze!
## Loading required package: sciplot
library(Rmisc)
library(sciplot)
library(plyr)
library(stringr) # for extracting numeric values
require(segmented) # for segmented line fitsß
## Loading required package: segmented
## Loading required package: MASS
## Loading required package: nlme
graz.data <- read.csv("/Users/hollyvm/GoogleSync/StudentWork/AGCamara/Grazing_all_Och_MH_MD - ExptData_v2.csv")
#graz.data<- graz.data[graz.data$OldData!= "Y",]
# From Barbaglia et al:
Barbaglia.data <- read.csv("/Users/hollyvm/GoogleSync/StudentWork/AGCamara/BarbagliaData/OchTradeoffData.csv")
# Rainbow of colours for the eight strains (as in Barbaglia et al 2024)
# Colours from: https://github.com/BlakeRMills/MetBrewer - Hokusai 1 + a bit of Hokusai 2
col.a <- rgb(96/255,26/255,16/255)
col.b <- rgb(183/255,51/255,48/255)
col.c <- rgb(234/255,95/255,75/255)
col.d <- rgb(233/255,120/255,47/255)
col.e <- rgb(243/255,186/255,81/255)
col.f <- rgb(120/255,171/255,125/255)
col.g <- rgb(0/255,93/255,148/255)
col.h <- rgb(0/255,56/255,98/255)
col.j <- rgb(64/255, 0/255, 128/255)
# Define the colors
col.a <- rgb(96/255,26/255,16/255)
col.b <- rgb(183/255,51/255,48/255)
col.c <- rgb(234/255,95/255,75/255)
col.d <- rgb(233/255,120/255,47/255)
col.e <- rgb(243/255,186/255,81/255)
col.f <- rgb(120/255,171/255,125/255)
col.g <- rgb(0/255,93/255,148/255)
col.h <- rgb(0/255,56/255,98/255)
col.j <- rgb(64/255, 0/255, 128/255)
# Rainbow of colours with reduced opacity for backgrounds of plots
alpha.bg <- .5
col.a.bg <- rgb(96/255,26/255,16/255,alpha.bg)
col.b.bg <- rgb(183/255,51/255,48/255,alpha.bg)
col.c.bg <- rgb(234/255,95/255,75/255,alpha.bg)
col.d.bg <- rgb(233/255,120/255,47/255,alpha.bg)
col.e.bg <- rgb(243/255,186/255,81/255,alpha.bg)
col.f.bg <- rgb(120/255,171/255,125/255,alpha.bg)
col.g.bg <- rgb(0/255,93/255,148/255,alpha.bg)
col.h.bg <- rgb(0/255,56/255,98/255,alpha.bg)
col.j.bg <- rgb(64/255, 0/255, 128/255, alpha.bg)
# Assembling into a vector of colors
col.strains <- c(col.a,col.b,col.c,col.d,col.e,col.f,col.g,col.h, col.j)
col.strains.bg <- c(col.a.bg,col.b.bg,col.c.bg,col.d.bg,col.e.bg,col.f.bg,col.g.bg,col.h.bg,col.j.bg)
# "Mega" vector alternates saturated and unsaturated, for making barplots that show high and low bacterial treatments side by side
mega.col.strains <- rep(NaN,length(col.strains)*2)
mega.col.strains.white <- rep(NaN,length(col.strains)*2)
for(i in 1:length(col.strains)){
mega.col.strains[2*(i-1)+1] <- col.strains.bg[i]
mega.col.strains.white[2*(i-1)+1] <- 'white'
mega.col.strains[2*i] <- col.strains[i]
mega.col.strains.white[2*i] <- col.strains[i]
}
# Line types for temperatures
ltys <- c(1,2)
# Shapes for replicates
rep.pchs <- c(1,4,5)
xenic.strains <- c(584,590,1148,1150,1391,1392,1393,2951)
rice.strains <- c('584R','590R','1148R','1150R','1391R',"1392R",'1393R','2951R')
# Make a list of all the different prey concentrations
preyconc <- levels(as.factor(graz.data$TargetOchConc))
preyconc
## [1] "0" "10000" "15000" "20000" "25000" "40000" "50000" "60000"
## [9] "70000" "75000" "90000" "95000" "100000" "106000" "125000"
# Dates for plotting
days <- c(1:4)
days
## [1] 1 2 3 4
# Select the columns to use for plotting
preycolumns <- c(12:15)+1
head(graz.data[,preycolumns]) # see how this subsets the data to just the prey columns?
## D1OchPop D2OchPop D3OchPop D4OchPop
## 1 741 228 969 4883
## 2 20511 15726 24557 33501
## 3 49282 31967 57264 63883
## 4 76973 48726 77484 87685
## 5 107225 76510 89041 85538
## 6 125564 105118 100437 102322
predcolumns <- c(17:20)+1
head(graz.data[,predcolumns])
## D1PredPop D2PredPop D3PredPop D4PredPop
## 1 0 0 0 0
## 2 0 0 0 0
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 0 0 0
## 6 0 0 0 0
# Set up colors for each prey concentration; using VanGogh3 palette from MetBrewer https://github.com/BlakeRMills/MetBrewer
ltgreen <- rgb(156/255,193/255,132/255)
dkgreen <- rgb(31/255,90/255,37/255)
ctrlcol <- 'gray80'
colfunc <- colorRampPalette(c(ltgreen,dkgreen))
col2use <- colfunc(length(preyconc))
col2use
## [1] "#9CC184" "#93B97D" "#8AB276" "#81AA6F" "#78A368" "#6F9C62" "#66945B"
## [8] "#5D8D54" "#54864D" "#4B7E46" "#427740" "#397039" "#306832" "#27612B"
## [15] "#1F5A25"
# Set up line types for different predation treatments
predtreats <- levels(as.factor(graz.data$Pred))
predtreats
## [1] "Control" "O.marina"
ltys <- c(1:length(predtreats))
# Set up shapes for different replicates
pchs <- c(21:25)
reps <- unique(graz.data$Rep)
# Create a list of unique experiments
strains <- sort(unique(graz.data$Strain))
Calculating growth rates
# Add new columns to our data frame that we will use for prey and predator growth rates and y-intercepts. We'll consider both 48 and 72hr windows for the prey.
graz.data$prey.mu.48 <- NaN
graz.data$prey.yint.48 <- NaN
graz.data$prey.R2.48 <- NaN
graz.data$prey.mu.72 <- NaN
graz.data$prey.yint.72 <- NaN
graz.data$prey.R2.72 <- NaN
graz.data$pred.mu.72 <- NaN
graz.data$pred.yint.72 <- NaN
graz.data$pred.R2.72 <- NaN
# Fit the growth curves line-by-line, using lm()
for( i in 1:dim(graz.data)[1] ){ # for each row in the data table
if(graz.data[i,]$AssayTemp %in% c(24)){ # excluding 30*C data
if(graz.data[i,]$TargetOchConc!=0){ # Only fit prey if we expect positive densities
# Fit all prey
lm1 <- lm(log(t(graz.data[i,preycolumns])+1)~days)
graz.data$prey.mu.72[i] <- summary(lm1)$coefficients[2,1]
graz.data$prey.yint.72[i] <- summary(lm1)$coefficients[1,1]
graz.data$prey.R2.72[i] <- summary(lm1)$r.squared
lm1 <- lm(log(t(graz.data[i,preycolumns[1:3]])+1)~days[1:3])
graz.data$prey.mu.48[i] <- summary(lm1)$coefficients[2,1]
graz.data$prey.yint.48[i] <- summary(lm1)$coefficients[1,1]
graz.data$prey.R2.48[i] <- summary(lm1)$r.squared
}
if( graz.data[i,]$Pred!='Control' ){ # Only fit predators if we expect non-zero densities
# Fit predator
if(!is.na(graz.data[i,]$D1PredPop)){
lm1 <- lm(log(t(graz.data[i,predcolumns]))~days)
graz.data$pred.mu.72[i] <- summary(lm1)$coefficients[2,1]
graz.data$pred.yint.72[i] <- summary(lm1)$coefficients[1,1]
graz.data$pred.R2.72[i] <- summary(lm1)$r.squared
}}
}
}
Let’s look at the differences in our 48 vs 72hr curve fits and check our work by adding our growth curves to our plots of population dynamics. This code creates one set of plots for every Ochromonas strain type; for strains for which multiple experiments were run, the data are broken down by experiment.
#par(mar=c(4,5,1,1),mfrow=c(4,2))
for(strainctr in 1:length(strains)){
#if(rem(strainctr,2)==1){
par(mar=c(4,5,1,1),mfrow=c(3,2))
# }
StrainChoice <- strains[strainctr]
AssayLightChoice <- 100
AssayTempChoice <- 24
# Subset the dataframe to the focal Ochromonas strain
dat5 <- graz.data[graz.data$Strain==StrainChoice & graz.data$AssayLight==AssayLightChoice&graz.data$AssayTemp==AssayTempChoice,]
if(strainctr==3){ dat5 <- dat5[dat5$OldData!='Y',]}
for(datectr in 1:length(unique(dat5$D1Date))){
dat1 <- dat5[dat5$D1Date==unique(dat5$D1Date)[datectr],]
# Create a vector of prey concentrations for this Ochromonas strain
preyconc <- levels(as.factor(dat1$TargetOchConc))
colfunc <- colorRampPalette(c('gray90',col.strains[ceiling(strainctr/2)]))
col2use <- colfunc(length(preyconc))
dat2 <- dat1[dat1$TargetOchConc!=0,]
plot(days,dat2[1,preycolumns]+1,las=1,xlab='Time (days)',ylab='',log='y',ylim=c((min(dat2[,preycolumns]+1,na.rm=T)),(max(dat2[,preycolumns]+1,na.rm=T))),type='n',main=paste('CCMP',StrainChoice, ', ', dat2$D1Date[1]))
mtext('Population size (cells/mL)',side=2,line=4)
for(i in 2:length(preyconc)){
subdat <- dat2[dat2$TargetOchConc==preyconc[i],]
for(j in 1:dim(subdat)[1]){
points(days,subdat[j,preycolumns],col=col2use[i],pch=21,bg=col2use[i])
if(subdat[j,]$Pred=='Control'){
# 48hr growth curve
xcoords <- seq(from = min(days), to = max(days[3]), length.out = 100)
ycoords <- exp(subdat[j,]$prey.yint.48)*exp(subdat[j,]$prey.mu.48*xcoords)
lines(xcoords,ycoords,col=col2use[i])
lines(xcoords,ycoords,col='red')
# 72hr growth curve
xcoords <- seq(from = min(days), to = max(days), length.out = 100)
ycoords <- exp(subdat[j,]$prey.yint.72)*exp(subdat[j,]$prey.mu.72*xcoords)
lines(xcoords,ycoords,col=col2use[i])
} else{
# 48hr growth curve
xcoords <- seq(from = min(days), to = max(days[3]), length.out = 100)
ycoords <- exp(subdat[j,]$prey.yint.48)*exp(subdat[j,]$prey.mu.48*xcoords)
lines(xcoords,ycoords,col=col2use[i],lty=2)
lines(xcoords,ycoords,col='red',lty=2)
# 72hr growth curve
xcoords <- seq(from = min(days), to = max(days), length.out = 100)
ycoords <- exp(subdat[j,]$prey.yint.72)*exp(subdat[j,]$prey.mu.72*xcoords)
lines(xcoords,ycoords,col=col2use[i],lty=2)
}
}
}
dat2 <- dat1[dat1$Pred!='Control',]
plot(days,dat2[1,predcolumns],las=1,xlab='Time (days)',ylab='Population size (cells/mL)',log='y',ylim=c((min(dat2[,predcolumns],na.rm=T)),(max(dat2[,predcolumns],na.rm=T))),type='n',main='Predator')
#mtext('Population size (cells/mL)',side=2,line=4)
for(i in 1:length(preyconc)){
subdat <- dat2[dat2$TargetOchConc==preyconc[i],]
for(j in 1:dim(subdat)[1]){
points(days,subdat[j,predcolumns],col=col2use[i],pch=21,bg=col2use[i])
xcoords <- seq(from = min(days), to = max(days), length.out = 100)
ycoords <- exp(subdat[j,]$pred.yint.72)*exp(subdat[j,]$pred.mu.72*xcoords)
lines(xcoords,ycoords,col=col2use[i])
}
}
}
}
Our goal is to use the growth rates of the different organisms to determine the grazing rate. A core part of this is determining differences in prey growth rates between experimental and control wells. We should see that, in the presence of O. marina, our prey grow slower.
The first thing that we need to do is match up our control (predator-free) growth rates with our predator ones.
graz.data$Strain.D1Date.TargetOchConc.Rep <- paste(graz.data$Strain,graz.data$D1Date,graz.data$TargetOchConc,graz.data$Rep,sep='.')
graz.data$ctrl.prey.mu.72 <- graz.data[graz.data$Pred=='Control',]$prey.mu.72[match(graz.data$Strain.D1Date.TargetOchConc.Rep,graz.data[graz.data$Pred=='Control',]$Strain.D1Date.TargetOchConc.Rep)] # the "match" function is a really powerful tool for linking data. Here, I'm using the "target prey" variable to make sure I match the correct average growth rates to each row in my data frame.
head(graz.data) # we can inspect the data frame to see if this worked. notice how the growth rates line up?
## OldData OldDat2 Strain Rep AssayTemp AssayLight TargetOchConc Pred D1Date
## 1 584 A 24 100 0 Control 1/30/24
## 2 584 A 24 100 25000 Control 1/30/24
## 3 584 A 24 100 50000 Control 1/30/24
## 4 584 A 24 100 75000 Control 1/30/24
## 5 584 A 24 100 100000 Control 1/30/24
## 6 584 A 24 100 125000 Control 1/30/24
## VolOch VolK X D1OchPop D2OchPop D3OchPop D4OchPop X.1 D1PredPop D2PredPop
## 1 NA NA NA 741 228 969 4883 NA 0 0
## 2 NA NA NA 20511 15726 24557 33501 NA 0 0
## 3 NA NA NA 49282 31967 57264 63883 NA 0 0
## 4 NA NA NA 76973 48726 77484 87685 NA 0 0
## 5 NA NA NA 107225 76510 89041 85538 NA 0 0
## 6 NA NA NA 125564 105118 100437 102322 NA 0 0
## D3PredPop D4PredPop X.2 Scope Run prey.mu.48 prey.yint.48 prey.R2.48
## 1 0 0 NA MD 1 NaN NaN NaN
## 2 0 0 NA MD 1 0.090013793 9.720203 0.1611998785
## 3 0 0 NA MD 1 0.075055211 10.560980 0.0614850373
## 4 0 0 NA MD 1 0.003308333 11.094400 0.0001547864
## 5 0 0 NA MD 1 -0.092915294 11.594080 0.3021313158
## 6 0 0 NA MD 1 -0.111641466 11.830191 0.8954060271
## prey.mu.72 prey.yint.72 prey.R2.72 pred.mu.72 pred.yint.72 pred.R2.72
## 1 NaN NaN NaN NaN NaN NaN
## 2 0.19174439 9.550652 0.6074022 NaN NaN NaN
## 3 0.13614234 10.459168 0.3344573 NaN NaN NaN
## 4 0.08547347 10.957459 0.1822492 NaN NaN NaN
## 5 -0.05262254 11.526925 0.2341172 NaN NaN NaN
## 6 -0.06596196 11.754058 0.6879644 NaN NaN NaN
## Strain.D1Date.TargetOchConc.Rep ctrl.prey.mu.72
## 1 584.1/30/24.0.A NaN
## 2 584.1/30/24.25000.A 0.19174439
## 3 584.1/30/24.50000.A 0.13614234
## 4 584.1/30/24.75000.A 0.08547347
## 5 584.1/30/24.100000.A -0.05262254
## 6 584.1/30/24.125000.A -0.06596196
However, because the control data are rather noisy for these experiments, we’re also going to use an alternate approach in which we use the mean growth rates from the control prey treatments and the mean growth rates for the predators. This helps to eliminate the rep-by-rep noise and give more accurate estimates of grazing rates.
graz.data$Strain.D1Date.TargetOchConc <- paste(graz.data$Strain,graz.data$D1Date,graz.data$TargetOchConc,sep='.')
ctrl.mu.summ <- summarySE(data=graz.data[graz.data$Pred=='Control',],measurevar='prey.mu.72',groupvars = 'Strain.D1Date.TargetOchConc')
graz.data$ctrl.prey.mu.72.mean <- ctrl.mu.summ$prey.mu.72[match(graz.data$Strain.D1Date.TargetOchConc,ctrl.mu.summ$Strain.D1Date.TargetOchConc)] # the "match" function is a really powerful tool for linking data. Here, I'm using the "target prey" variable to make sure I match the correct average growth rates to each row in my data frame.
head(graz.data) # we can inspect the data frame to see if this worked. notice how the growth rates line up?
## OldData OldDat2 Strain Rep AssayTemp AssayLight TargetOchConc Pred D1Date
## 1 584 A 24 100 0 Control 1/30/24
## 2 584 A 24 100 25000 Control 1/30/24
## 3 584 A 24 100 50000 Control 1/30/24
## 4 584 A 24 100 75000 Control 1/30/24
## 5 584 A 24 100 100000 Control 1/30/24
## 6 584 A 24 100 125000 Control 1/30/24
## VolOch VolK X D1OchPop D2OchPop D3OchPop D4OchPop X.1 D1PredPop D2PredPop
## 1 NA NA NA 741 228 969 4883 NA 0 0
## 2 NA NA NA 20511 15726 24557 33501 NA 0 0
## 3 NA NA NA 49282 31967 57264 63883 NA 0 0
## 4 NA NA NA 76973 48726 77484 87685 NA 0 0
## 5 NA NA NA 107225 76510 89041 85538 NA 0 0
## 6 NA NA NA 125564 105118 100437 102322 NA 0 0
## D3PredPop D4PredPop X.2 Scope Run prey.mu.48 prey.yint.48 prey.R2.48
## 1 0 0 NA MD 1 NaN NaN NaN
## 2 0 0 NA MD 1 0.090013793 9.720203 0.1611998785
## 3 0 0 NA MD 1 0.075055211 10.560980 0.0614850373
## 4 0 0 NA MD 1 0.003308333 11.094400 0.0001547864
## 5 0 0 NA MD 1 -0.092915294 11.594080 0.3021313158
## 6 0 0 NA MD 1 -0.111641466 11.830191 0.8954060271
## prey.mu.72 prey.yint.72 prey.R2.72 pred.mu.72 pred.yint.72 pred.R2.72
## 1 NaN NaN NaN NaN NaN NaN
## 2 0.19174439 9.550652 0.6074022 NaN NaN NaN
## 3 0.13614234 10.459168 0.3344573 NaN NaN NaN
## 4 0.08547347 10.957459 0.1822492 NaN NaN NaN
## 5 -0.05262254 11.526925 0.2341172 NaN NaN NaN
## 6 -0.06596196 11.754058 0.6879644 NaN NaN NaN
## Strain.D1Date.TargetOchConc.Rep ctrl.prey.mu.72 Strain.D1Date.TargetOchConc
## 1 584.1/30/24.0.A NaN 584.1/30/24.0
## 2 584.1/30/24.25000.A 0.19174439 584.1/30/24.25000
## 3 584.1/30/24.50000.A 0.13614234 584.1/30/24.50000
## 4 584.1/30/24.75000.A 0.08547347 584.1/30/24.75000
## 5 584.1/30/24.100000.A -0.05262254 584.1/30/24.100000
## 6 584.1/30/24.125000.A -0.06596196 584.1/30/24.125000
## ctrl.prey.mu.72.mean
## 1 NaN
## 2 0.17969735
## 3 0.15188058
## 4 0.05697058
## 5 0.02953099
## 6 -0.01727983
From here, we can compute our grazing rates as the difference between the control and +predator growth rates.
# Rep-by-rep
graz.data$prey.g.72 <- graz.data$ctrl.prey.mu.72-graz.data$prey.mu.72
graz.data$prey.g.72 <- graz.data$ctrl.prey.mu.72-graz.data$prey.mu.72
# Mean
graz.data$prey.g.72.mean <- graz.data$ctrl.prey.mu.72.mean-graz.data$prey.mu.72
graz.data$prey.g.72.mean <- graz.data$ctrl.prey.mu.72.mean-graz.data$prey.mu.72
Computing the per-predator ingestion rate is a little trickier. Ask someone to show you the calculus on the whiteboard, but essentially what we are doing is using the growth rates of both predators and prey to compute mean populations of both, and then computing the per-predator ingestion rate from there. This method follows Jeong & Latz (1994) (https://www.int-res.com/articles/meps/106/m106p173.pdf)
graz.data$prey.popn.24 <- graz.data$D1OchPop*(exp(graz.data$prey.mu.48)-1)/(graz.data$prey.mu.48) # mean prey population
graz.data[graz.data$TargetOchConc==0,]$prey.popn.24 <- 0
graz.data$pred.popn.24 <- graz.data$D1PredPop*(exp(graz.data$pred.mu.72)-1)/(graz.data$pred.mu.72) # mean predator population
graz.data$prey.i.24 <- graz.data$prey.g.72*graz.data$prey.popn.24/graz.data$pred.popn.24 # ingestion = mean prey * grazing rate / mean predator
graz.data$prey.i.24.mean <- graz.data$prey.g.72.mean*graz.data$prey.popn.24/graz.data$pred.popn.24 # ingestion = mean prey * grazing rate / mean predator
# Finally we'll manually fill in grazing and ingestion rates of 0 for 0 prey treatments.
graz.data[graz.data$TargetOchConc==0&graz.data$Pred!='Control',]$prey.g.72 <- 0
graz.data[graz.data$TargetOchConc==0&graz.data$Pred!='Control',]$prey.i.24 <- 0
graz.data[graz.data$TargetOchConc==0&graz.data$Pred!='Control',]$prey.g.72.mean <- 0
graz.data[graz.data$TargetOchConc==0&graz.data$Pred!='Control',]$prey.i.24.mean <- 0
Let’s look at our results! A scatterplot lets us see our results as a function of the mean prey concentration. We can even work on fitting our data with a Holling Type II (saturating) functional response. Here, I’m also choosing to fit a line to the most linear part of the data to estimate the predator’s attack rate.
par(mar=c(2,2.5,1,1),mfrow=c(2,4))
linear.fit.maxima <- c(70000, 50000, 40000, 25000, 155000,90000, 60000, 80000)
saturating.fit.maxima <- c(100000, 100000, 90000, 135000, 155000,250000, 170000, 110000)
attack.ests <- as.data.frame(xenic.strains)
colnames(attack.ests) <-'Strain'
attack.ests$a.linear <- NaN
attack.ests$a.linear.se <- NaN
attack.ests$a.sat <- NaN
attack.ests$a.sat.se <- NaN
attack.ests$h.sat <- NaN
attack.ests$h.sat.se <- NaN
for(i in 1:length(xenic.strains)){
subdat1 <- rbind(graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
# Make the plot
plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.popn.24,na.rm=T)), ylim=c(0,max(subdat1$prey.i.24,na.rm=T)),las=1)
# Add labeling
text(x = 0-.06*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$prey.i.24,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)
if(i == 1){
legend(x='topright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2)
}
# Add xenic data
subdat.x <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
points(subdat.x$prey.popn.24,subdat.x$prey.i.24,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
# Add rice data
subdat.r <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
points(subdat.r$prey.popn.24,subdat.r$prey.i.24, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
# Add linear model for all data
subdat1.trim <- subdat1[subdat1$prey.popn.24 < linear.fit.maxima[i],]
lm1 <- nls(prey.i.24 ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
xset <- seq(from = 0, to = linear.fit.maxima[i])
yset <- summary(lm1)$parameters[1] * xset
lines(xset,yset,lwd=1.5,col=col.strains[i])
#abline(lm1)
attack.ests$a.linear[i] <- summary(lm1)$parameters[1]
attack.ests$a.linear.se[i] <- summary(lm1)$parameters[2]
# Add saturating model for all data
subdat1.trim <- subdat1[subdat1$prey.popn.24 < saturating.fit.maxima[i],]
a.est <- summary(lm(subdat1.trim[1:7,]$prey.i.24~subdat1.trim[1:7,]$prey.popn.24))$coefficients[2,1] # estimate initial slope
h.est <- 1/50 # estimate handling time
fit1 <- nls(prey.i.24 ~ a*prey.popn.24/(1+a*h*prey.popn.24),data=subdat1.trim,start = list(a = a.est, h = h.est))
a.fit <- summary(fit1)$parameters[1,1]
h.fit <- summary(fit1)$parameters[2,1]
xcoords <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T),length.out = 100)
ycoords <- a.fit*xcoords/(1+a.fit*h.fit*xcoords)
lines(xcoords,ycoords,lwd=2,col='black',lty=2)
attack.ests$a.sat[i] <- summary(fit1)$parameters[1,1]
attack.ests$a.sat.se[i] <- summary(fit1)$parameters[1,2]
attack.ests$h.sat[i] <- summary(fit1)$parameters[2,1]
attack.ests$h.sat.se[i] <- summary(fit1)$parameters[2,2]
}
Fitting linear models with cutoffs
par(mar=c(2,2.5,1,1),mfrow=c(2,4))
linear.fit.maxima <- c(70000, 50000, 40000, 25000, 155000,90000, 60000, 80000)
linear.fit.maxima.x <- c(150000, 100000, 50000, 25000, 155000,120000, 600000, 120000)
linear.fit.maxima.r <- c(70000, 100000, 40000, 150000, 155000,250000, 200000, 120000)
saturating.fit.maxima <- c(100000, 100000, 90000, 135000, 155000,250000, 170000, 110000)
attack.ests <- as.data.frame(xenic.strains)
colnames(attack.ests) <-'Strain'
attack.ests$a.linear <- NaN
attack.ests$a.linear.se <- NaN
attack.ests$a.x.linear <- NaN
attack.ests$a.x.linear.se <- NaN
attack.ests$a.r.linear <- NaN
attack.ests$a.r.linear.se <- NaN
attack.ests$a.sat <- NaN
attack.ests$a.sat.se <- NaN
attack.ests$h.sat <- NaN
attack.ests$h.sat.se <- NaN
for(i in 1:length(xenic.strains)){
subdat1 <- rbind(graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
# Make the plot
plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.popn.24,na.rm=T)), ylim=c(0,max(subdat1$prey.i.24,na.rm=T)),las=1)
# Add labeling
text(x = 0-.06*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$prey.i.24,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)
if(i == 1){
legend(x='topright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2)
}
# Add xenic data
subdat.x <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
points(subdat.x$prey.popn.24,subdat.x$prey.i.24,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
# Add rice data
subdat.r <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
points(subdat.r$prey.popn.24,subdat.r$prey.i.24, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
# Add linear model for all data
subdat1.trim <- subdat1[subdat1$prey.popn.24 < linear.fit.maxima[i],]
lm1 <- nls(prey.i.24 ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
xset <- seq(from = 0, to = linear.fit.maxima[i])
yset <- summary(lm1)$parameters[1] * xset
#lines(xset,yset,lwd=1.5,col=col.strains[i])
#abline(lm1)
attack.ests$a.linear[i] <- summary(lm1)$parameters[1]
attack.ests$a.linear.se[i] <- summary(lm1)$parameters[2]
# Add linear model for xenic data
subdat1.trim <- subdat1[subdat1$Strain==xenic.strains[i] ,]
subdat1.trim <- subdat1[subdat1$Strain==xenic.strains[i] & subdat1$prey.popn.24 < linear.fit.maxima.x[i],]
lm1 <- nls(prey.i.24 ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
xset <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T),length.out=100)
yset <- summary(lm1)$parameters[1] * xset
lines(xset,yset,lwd=1.5,col=col.strains[i],lty=2)
attack.ests$a.x.linear[i] <- summary(lm1)$parameters[1]
attack.ests$a.x.linear.se[i] <- summary(lm1)$parameters[2]
# Add linear model for rice data
subdat1.trim <- subdat1[subdat1$Strain==rice.strains[i],]
subdat1.trim <- subdat1[subdat1$Strain==rice.strains[i] & subdat1$prey.popn.24 < linear.fit.maxima.r[i],]
lm1 <- nls(prey.i.24 ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
xset <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T))
yset <- summary(lm1)$parameters[1] * xset
lines(xset,yset,lwd=1.5,col=col.strains[i],lty=1)
attack.ests$a.r.linear[i] <- summary(lm1)$parameters[1]
attack.ests$a.r.linear.se[i] <- summary(lm1)$parameters[2]
# Add saturating model for all data
subdat1.trim <- subdat1[subdat1$prey.popn.24 < saturating.fit.maxima[i],]
a.est <- summary(lm(subdat1.trim[1:7,]$prey.i.24~subdat1.trim[1:7,]$prey.popn.24))$coefficients[2,1] # estimate initial slope
h.est <- 1/50 # estimate handling time
fit1 <- nls(prey.i.24 ~ a*prey.popn.24/(1+a*h*prey.popn.24),data=subdat1.trim,start = list(a = a.est, h = h.est))
a.fit <- summary(fit1)$parameters[1,1]
h.fit <- summary(fit1)$parameters[2,1]
xcoords <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T),length.out = 100)
ycoords <- a.fit*xcoords/(1+a.fit*h.fit*xcoords)
#lines(xcoords,ycoords,lwd=2,col='black',lty=2)
attack.ests$a.sat[i] <- summary(fit1)$parameters[1,1]
attack.ests$a.sat.se[i] <- summary(fit1)$parameters[1,2]
attack.ests$h.sat[i] <- summary(fit1)$parameters[2,1]
attack.ests$h.sat.se[i] <- summary(fit1)$parameters[2,2]
}
Adjust with breakpoint analysis to make cutoffs more statistically robust
par(mar=c(2,2.5,.7,.5))
rowsperpanel <- 5
colsperpanel <- 5
lay.mat <- matrix(data=NA,nrow=2*rowsperpanel+1,ncol=4*colsperpanel+1)
lay.mat[,1] <- 9
lay.mat[2*rowsperpanel+1,] <- 10
for(i in 1:2){
for(j in 1:4){
lay.mat[(1+(i-1)*rowsperpanel):(i*rowsperpanel),(2+(j-1)*colsperpanel):(1+j*colsperpanel)] <- 4*(i-1) + j
}
}
layout(lay.mat)
# Create a data frame to hold the results
attack.ests <- as.data.frame(xenic.strains)
colnames(attack.ests) <-'Strain'
attack.ests$a.linear <- NaN
attack.ests$a.linear.se <- NaN
attack.ests$a.x.linear <- NaN
attack.ests$a.x.linear.se <- NaN
attack.ests$a.r.linear <- NaN
attack.ests$a.r.linear.se <- NaN
attack.ests$a.sat <- NaN
attack.ests$a.sat.se <- NaN
attack.ests$h.sat <- NaN
attack.ests$h.sat.se <- NaN
for(i in 1:length(xenic.strains)){
subdat1 <- rbind(graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
# Make the plot
plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.popn.24,na.rm=T)), ylim=c(0,max(subdat1$prey.i.24,na.rm=T)),las=1)
# Add labeling
text(x = 0-.06*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$prey.i.24,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)
if(i == 8){
legend(x='topright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2)
}
if(i == 1){
mtext(expression(paste('Ingestion rate (',italic('Ochromonas'),' · ',italic('O. marina')^-1,' · ',d^-1,')')),side=2,cex=.8,at=-6,line=2.5)
}
if(i == 6){
mtext(expression(paste('Prey availability (',italic('Ochromonas'),' · ',mL^-1, ')')),side=1,line=2.7,cex=.8,at=275000)
}
# Add xenic data
subdat.x <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
points(subdat.x$prey.popn.24,subdat.x$prey.i.24,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
# Add rice data
subdat.r <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
points(subdat.r$prey.popn.24,subdat.r$prey.i.24, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
# Add linear model for all data
## Test 3 types of model fits
lm1 <- lm(prey.i.24 ~ prey.popn.24+0, data=subdat1) # +0 forces y-intercept through 0
o1 <- segmented(lm1)
o2 <- segmented(lm1,npsi=2)
## Select amongst them based on smallest p-value
seg.pvals <- c(summary(lm1)$coefficients[1,4],summary(o1)$coefficients[1,4],summary(o2)$coefficients[1,4])
if(which(seg.pvals==min(seg.pvals))[1]==1){ # If the standard linear model is the best
# Plot the best-fit model; skipping here
# Record the data
attack.ests$a.linear[i] <- summary(lm1)$coefficients[1,1]
attack.ests$a.linear.se[i] <- summary(lm1)$coefficients[1,2]
}
if(which(seg.pvals==min(seg.pvals))[1]==2){ # If the 1 breakpoint model is the best
# Plot the best-fit model; skipping here
# Record the data
attack.ests$a.linear[i] <- summary(o1)$coefficients[1,1]
attack.ests$a.linear.se[i] <- summary(o1)$coefficients[1,2]
}
if(which(seg.pvals==min(seg.pvals))[1]==3){ # If the 2 breakpoint model is the best
# Plot the best-fit model; skipping here
# Record the data
attack.ests$a.linear[i] <- summary(o2)$coefficients[1,1]
attack.ests$a.linear.se[i] <- summary(o2)$coefficients[1,2]
}
# Add linear model for xenic data
## Test 3 types of model fits
lm1 <- lm(prey.i.24 ~ prey.popn.24+0, data=subdat1[subdat1$Strain==xenic.strains[i] ,])
o1 <- segmented(lm1)
o2 <- segmented(lm1,npsi=2)
## Select amongst them based on smallest p-value
seg.pvals <- c(summary(lm1)$coefficients[1,4],summary(o1)$coefficients[1,4],summary(o2)$coefficients[1,4])
if(which(seg.pvals==min(seg.pvals))[1]==1){ # If the standard linear model is the best
# Plot the best-fit model
xcoords <- seq(from = 0, to = max(subdat1[subdat1$Strain==xenic.strains[i],]$prey.popn.24,na.rm=T),length.out=100)
ycoords <- xcoords*summary(lm1)$coefficients[1,1]
lines(xcoords,ycoords,lwd=2,lty=2,col=col.strains[i])
# Record the data
attack.ests$a.x.linear[i] <- summary(lm1)$coefficients[1,1]
attack.ests$a.x.linear.se[i] <- summary(lm1)$coefficients[1,2]
}
if(which(seg.pvals==min(seg.pvals))[1]==2){ # If the 1 breakpoint model is the best
# Plot the best-fit model
plot.segmented(o1,add=TRUE,col=col.strains[i],lty=2,lwd=2)
# Record the data
attack.ests$a.x.linear[i] <- summary(o1)$coefficients[1,1]
attack.ests$a.x.linear.se[i] <- summary(o1)$coefficients[1,2]
}
if(which(seg.pvals==min(seg.pvals))[1]==3){ # If the 2 breakpoint model is the best
# Plot the best-fit model
plot.segmented(o2,add=TRUE,col=col.strains[i],lty=2,lwd=2)
# Record the data
attack.ests$a.x.linear[i] <- summary(o2)$coefficients[1,1]
attack.ests$a.x.linear.se[i] <- summary(o2)$coefficients[1,2]
}
# Add linear model for rice data
## Test 3 types of model fits
lm1 <- lm(prey.i.24 ~ prey.popn.24+0, data=subdat1[subdat1$Strain==rice.strains[i] ,])
o1 <- segmented(lm1)
if(i == 6){
o2 <- o1
} else{
o2 <- segmented(lm1,npsi=2)}
## Select amongst them based on smallest p-value
seg.pvals <- c(summary(lm1)$coefficients[1,4],summary(o1)$coefficients[1,4],summary(o2)$coefficients[1,4])
if(which(seg.pvals==min(seg.pvals))[1]==1){ # If the standard linear model is the best
# Plot the best-fit model
xcoords <- seq(from = 0, to = max(subdat1[subdat1$Strain==rice.strains[i],]$prey.popn.24,na.rm=T),length.out=100)
ycoords <- xcoords*summary(lm1)$coefficients[1,1]
lines(xcoords,ycoords,lwd=2,lty=1,col=col.strains[i])
# Record the data
attack.ests$a.r.linear[i] <- summary(lm1)$coefficients[1,1]
attack.ests$a.r.linear.se[i] <- summary(lm1)$coefficients[1,2]
}
if(which(seg.pvals==min(seg.pvals))[1]==2){ # If the 1 breakpoint model is the best
# Plot the best-fit model
plot.segmented(o1,add=TRUE,col=col.strains[i],lty=1,lwd=2)
# Record the data
attack.ests$a.r.linear[i] <- summary(o1)$coefficients[1,1]
attack.ests$a.r.linear.se[i] <- summary(o1)$coefficients[1,2]
}
if(which(seg.pvals==min(seg.pvals))[1]==3){ # If the 2 breakpoint model is the best
# Plot the best-fit model
plot.segmented(o2,add=TRUE,col=col.strains[i],lty=1,lwd=2)
# Record the data
attack.ests$a.r.linear[i] <- summary(o2)$coefficients[1,1]
attack.ests$a.r.linear.se[i] <- summary(o2)$coefficients[1,2]
}
}
Including just the first line segment
par(mar=c(2,2.5,.7,.5))
rowsperpanel <- 5
colsperpanel <- 5
lay.mat <- matrix(data=NA,nrow=2*rowsperpanel+1,ncol=4*colsperpanel+1)
lay.mat[,1] <- 9
lay.mat[2*rowsperpanel+1,] <- 10
for(i in 1:2){
for(j in 1:4){
lay.mat[(1+(i-1)*rowsperpanel):(i*rowsperpanel),(2+(j-1)*colsperpanel):(1+j*colsperpanel)] <- 4*(i-1) + j
}
}
layout(lay.mat)
# Create a data frame to hold the results
attack.ests <- as.data.frame(xenic.strains)
colnames(attack.ests) <-'Strain'
attack.ests$a.linear <- NaN
attack.ests$a.linear.se <- NaN
attack.ests$a.x.linear <- NaN
attack.ests$a.x.linear.se <- NaN
attack.ests$a.r.linear <- NaN
attack.ests$a.r.linear.se <- NaN
attack.ests$a.sat <- NaN
attack.ests$a.sat.se <- NaN
attack.ests$h.sat <- NaN
attack.ests$h.sat.se <- NaN
for(i in 1:length(xenic.strains)){
subdat1 <- rbind(graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
# Make the plot
plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.popn.24,na.rm=T)), ylim=c(0,max(subdat1$prey.i.24,na.rm=T)),las=1)
# Add labeling
text(x = 0-.06*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$prey.i.24,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)
if(i == 1){
legend(x='topright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2)
mtext(expression(paste('Ingestion rate (',italic('Ochromonas'),' · ',italic('O. marina')^-1,' · ',d^-1,')')),side=2,cex=.8,at=-6,line=2.5)
}
if(i == 6){
mtext(expression(paste('Prey availability (',italic('Ochromonas'),' · ',mL^-1, ')')),side=1,line=2.7,cex=.8,at=275000)
}
# Add xenic data
subdat.x <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
points(subdat.x$prey.popn.24,subdat.x$prey.i.24,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
# Add rice data
subdat.r <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
points(subdat.r$prey.popn.24,subdat.r$prey.i.24, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
# Add linear model for all data
## Test 3 types of model fits
lm1 <- lm(prey.i.24 ~ prey.popn.24+0, data=subdat1) # +0 forces y-intercept through 0
o1 <- segmented(lm1)
o2 <- segmented(lm1,npsi=2)
## Select amongst them based on smallest p-value
seg.pvals <- c(summary(lm1)$coefficients[1,4],summary(o1)$coefficients[1,4],summary(o2)$coefficients[1,4])
if(which(seg.pvals==min(seg.pvals))[1]==1){ # If the standard linear model is the best
# Plot the best-fit model; skipping here
# Record the data
attack.ests$a.linear[i] <- summary(lm1)$coefficients[1,1]
attack.ests$a.linear.se[i] <- summary(lm1)$coefficients[1,2]
}
if(which(seg.pvals==min(seg.pvals))[1]==2){ # If the 1 breakpoint model is the best
# Plot the best-fit model; skipping here
# Record the data
attack.ests$a.linear[i] <- summary(o1)$coefficients[1,1]
attack.ests$a.linear.se[i] <- summary(o1)$coefficients[1,2]
}
if(which(seg.pvals==min(seg.pvals))[1]==3){ # If the 2 breakpoint model is the best
# Plot the best-fit model; skipping here
# Record the data
attack.ests$a.linear[i] <- summary(o2)$coefficients[1,1]
attack.ests$a.linear.se[i] <- summary(o2)$coefficients[1,2]
}
# Add linear model for xenic data
## Test 3 types of model fits
lm1 <- lm(prey.i.24 ~ prey.popn.24+0, data=subdat1[subdat1$Strain==xenic.strains[i] ,])
o1 <- segmented(lm1)
o2 <- segmented(lm1,npsi=2)
## Select amongst them based on smallest p-value
seg.pvals <- c(summary(lm1)$coefficients[1,4],summary(o1)$coefficients[1,4],summary(o2)$coefficients[1,4])
if(which(seg.pvals==min(seg.pvals))[1]==1){ # If the standard linear model is the best
# Plot the best-fit model
xcoords <- seq(from = 0, to = max(subdat1[subdat1$Strain==xenic.strains[i],]$prey.popn.24,na.rm=T),length.out=100)
ycoords <- xcoords*summary(lm1)$coefficients[1,1]
lines(xcoords,ycoords,lwd=2.5,lty=2,col=col.strains[i])
# Record the data
attack.ests$a.x.linear[i] <- summary(lm1)$coefficients[1,1]
attack.ests$a.x.linear.se[i] <- summary(lm1)$coefficients[1,2]
}
if(which(seg.pvals==min(seg.pvals))[1]==2){ # If the 1 breakpoint model is the best
# Plot the best-fit model
#plot.segmented(o1,add=TRUE,col=col.strains[i],lty=2,lwd=2)
xcoords <- seq(from = 0, to = o1$psi[1,2],length.out=100)
ycoords <- xcoords*summary(o1)$coefficients[1,1]
lines(xcoords,ycoords,lwd=2.5,lty=2,col=col.strains[i])
# Record the data
attack.ests$a.x.linear[i] <- summary(o1)$coefficients[1,1]
attack.ests$a.x.linear.se[i] <- summary(o1)$coefficients[1,2]
}
if(which(seg.pvals==min(seg.pvals))[1]==3){ # If the 2 breakpoint model is the best
# Plot the best-fit model
xcoords <- seq(from = 0, to = o2$psi[1,2],length.out=100)
ycoords <- xcoords*summary(o2)$coefficients[1,1]
lines(xcoords,ycoords,lwd=2.5,lty=2,col=col.strains[i])
# Record the data
attack.ests$a.x.linear[i] <- summary(o2)$coefficients[1,1]
attack.ests$a.x.linear.se[i] <- summary(o2)$coefficients[1,2]
}
# Add linear model for rice data
## Test 3 types of model fits
lm1 <- lm(prey.i.24 ~ prey.popn.24+0, data=subdat1[subdat1$Strain==rice.strains[i] ,])
o1 <- segmented(lm1)
if(i == 6){
o2 <- o1
} else{
o2 <- segmented(lm1,npsi=2)}
## Select amongst them based on smallest p-value
seg.pvals <- c(summary(lm1)$coefficients[1,4],summary(o1)$coefficients[1,4],summary(o2)$coefficients[1,4])
if(which(seg.pvals==min(seg.pvals))[1]==1){ # If the standard linear model is the best
# Plot the best-fit model
xcoords <- seq(from = 0, to = max(subdat1[subdat1$Strain==rice.strains[i],]$prey.popn.24,na.rm=T),length.out=100)
ycoords <- xcoords*summary(lm1)$coefficients[1,1]
lines(xcoords,ycoords,lwd=2,lty=1,col=col.strains[i])
# Record the data
attack.ests$a.r.linear[i] <- summary(lm1)$coefficients[1,1]
attack.ests$a.r.linear.se[i] <- summary(lm1)$coefficients[1,2]
}
if(which(seg.pvals==min(seg.pvals))[1]==2){ # If the 1 breakpoint model is the best
# Plot the best-fit model
#plot.segmented(o1,add=TRUE,col=col.strains[i],lty=2,lwd=2)
xcoords <- seq(from = 0, to = o1$psi[1,2],length.out=100)
ycoords <- xcoords*summary(o1)$coefficients[1,1]
lines(xcoords,ycoords,lwd=2,lty=1,col=col.strains[i])
# Record the data
attack.ests$a.r.linear[i] <- summary(o1)$coefficients[1,1]
attack.ests$a.r.linear.se[i] <- summary(o1)$coefficients[1,2]
}
if(which(seg.pvals==min(seg.pvals))[1]==3){ # If the 2 breakpoint model is the best
# Plot the best-fit model
xcoords <- seq(from = 0, to = o2$psi[1,2],length.out=100)
ycoords <- xcoords*summary(o2)$coefficients[1,1]
lines(xcoords,ycoords,lwd=2,lty=1,col=col.strains[i])
# Record the data
attack.ests$a.r.linear[i] <- summary(o2)$coefficients[1,1]
attack.ests$a.r.linear.se[i] <- summary(o2)$coefficients[1,2]
}
}
par(mar=c(4,4.2,1,1))
barplot(attack.ests$a.linear*1000,las=1, ylim=c(0,max(attack.ests$a.linear+attack.ests$a.linear.se))*1001,names=xenic.strains,xlab='Strain',ylab=expression(paste(italic('O. marina'),' clearance rate (μL ',d^-1,")")), border="black", col = col.strains.bg)
xcoords <- seq(from = 0.75, to = 9.1, length.out=8)
arrows(xcoords,attack.ests$a.linear*1000-attack.ests$a.linear.se*1000,xcoords, attack.ests$a.linear*1000+attack.ests$a.linear.se*1000, code = 3, angle = 90, length = 0.05, col= 'black')
subdat <- as.data.frame(cbind(attack.ests$Strain,attack.ests$a.r.linear,attack.ests$a.r.linear.se))
colnames(subdat) <- c('Strain','attack','se')
subdat$Bact <- 'R'
subdat2 <- as.data.frame(cbind(attack.ests$Strain,attack.ests$a.x.linear,attack.ests$a.x.linear.se))
colnames(subdat2) <- c('Strain','attack','se')
subdat2$Bact <- 'X'
subdat <- rbind(subdat2,subdat)
subdat <- subdat[order(subdat$Strain),]
barplot(subdat$attack*1000,las=1, ylim=c(0,max(subdat$attack+subdat$se))*1001,xlab='Strain',ylab=expression(paste(italic('O. marina'),' clearance rate (μL · ',d^-1,")")), border=rep(col.strains,each=2),lwd=1, col = mega.col.strains.white,space=c(0.7,.1))
xcoords1 <- seq(from = 1.2, to = 20.8, length.out=8)
xcoords2 <- seq(from = 2.3, to = 21.9, length.out = 8)
xcoords <- sort(c(xcoords1,xcoords2))
arrows(xcoords,subdat$attack*1000-subdat$se*1000,xcoords, subdat$attack*1000+subdat$se*1000, code = 3, angle = 90, length = 0.05, col= 'black')
axis(side=1, at = colMeans(rbind(xcoords1,xcoords2)), labels=xenic.strains, padj=-.5 )
legend(x='topright',inset=c(0,0),legend=c('Low Bacteria','High Bacteria'),pch=22,pt.cex=2,bty='n',pt.bg=c('white', 'gray10'), col=c('gray10','gray10'))
# Adding t-test results
barplot(subdat$attack*1000,las=1, ylim=c(0,max(subdat$attack+subdat$se))*1001,xlab='Strain',ylab=expression(paste(italic('O. marina'),' clearance rate (μL · ',d^-1,")")), border=rep(col.strains,each=2),lwd=1, col = mega.col.strains.white,space=c(0.7,.1))
xcoords1 <- seq(from = 1.2, to = 20.8, length.out=8)
xcoords2 <- seq(from = 2.3, to = 21.9, length.out = 8)
xcoords <- sort(c(xcoords1,xcoords2))
xcoords.mean <- colMeans(rbind(xcoords1,xcoords2))
arrows(xcoords,subdat$attack*1000-subdat$se*1000,xcoords, subdat$attack*1000+subdat$se*1000, code = 3, angle = 90, length = 0.05, col= 'black')
axis(side=1, at = xcoords.mean, labels=xenic.strains, padj=-.5 )
legend(x='topright',inset=c(0,0),legend=c('Low Bacteria','High Bacteria'),pch=22,pt.cex=2,bty='n',pt.bg=c('white', 'gray10'), col=c('gray10','gray10'))
pval.text <- rep('',8)
for(i in 1:length(xenic.strains)){
subdat2 <- subdat[subdat$Strain==xenic.strains[i],]
max.se <- max(subdat2$se)
delta.attack <- max(subdat2$attack) - min(subdat2$attack)
num.ses <- delta.attack / max.se
if(num.ses > 3){
pval.text[i] <- '*'
}
if(num.ses > 4){
pval.text[i] <- '**'
}
if(num.ses > 5){
pval.text[i] <- '***'
}
}
ycoords <- tapply(subdat$attack*1000,subdat$Strain,FUN='max') + .04
text(xcoords.mean,ycoords,pval.text,cex=.75)
Raw growth rates
par(mar=c(2,2.5,.7,.5))
rowsperpanel <- 5
colsperpanel <- 5
lay.mat <- matrix(data=NA,nrow=2*rowsperpanel+1,ncol=4*colsperpanel+1)
lay.mat[,1] <- 9
lay.mat[2*rowsperpanel+1,] <- 10
for(i in 1:2){
for(j in 1:4){
lay.mat[(1+(i-1)*rowsperpanel):(i*rowsperpanel),(2+(j-1)*colsperpanel):(1+j*colsperpanel)] <- 4*(i-1) + j
}
}
layout(lay.mat)
for(i in 1:length(xenic.strains)){
subdat1 <- rbind(graz.data[graz.data$TargetOchConc!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
# Make the plot
plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.popn.24,na.rm=T)), ylim=c(min(subdat1$pred.mu.72,na.rm=T),max(subdat1$pred.mu.72,na.rm=T)),las=1)
abline(h=0,lty=2)
# Add xenic data
subdat.x <- subdat1[subdat1$Strain==xenic.strains[i] ,]
points(subdat.x$prey.popn.24,subdat.x$pred.mu.72,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
# Add rice data
subdat.r <- subdat1[subdat1$Strain==rice.strains[i],]
points(subdat.r$prey.popn.24,subdat.r$pred.mu.72, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
# Add labeling
text(x = 0-.06*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$pred.mu.72,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)
if(i == 1){
legend(x='topright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2)
}
if(i == 1){
mtext(expression(paste(italic('O. marina'),' growth rate (',d^-1,')')),side=2,cex=.8,at=-.6,line=2.5)
}
if(i == 6){
mtext(expression(paste('Prey availability (',italic('Ochromonas'),' · ',mL^-1, ')')),side=1,line=2.7,cex=.8,at=275000)
}
}
Relativized growth rates
# First, calculate the mean growth rates of O. marina for each experiment
# Create a unique identifier for each experiment
graz.data$Strain.D1Date <- paste(graz.data$Strain,graz.data$D1Date,sep='.')
# Compute O. marina control (0 prey) growth rates for each experiment
ref.Omarina.growthrates <- summarySE(data=graz.data[graz.data$TargetOchConc==0 & !is.na(graz.data$pred.mu.72),],'pred.mu.72',groupvars='Strain.D1Date')
# Second, map the O. marina control growth rates back to their relevant experiment
graz.data$pred.mu.72.ctrl <- ref.Omarina.growthrates$pred.mu.72[match(graz.data$Strain.D1Date,ref.Omarina.growthrates$Strain.D1Date)]
# Third, compute the relativized growth rates
graz.data$pred.mu.72.rel <- graz.data$pred.mu.72-graz.data$pred.mu.72.ctrl
par(mar=c(2,2.5,.7,.5))
rowsperpanel <- 5
colsperpanel <- 5
lay.mat <- matrix(data=NA,nrow=2*rowsperpanel+1,ncol=4*colsperpanel+1)
lay.mat[,1] <- 9
lay.mat[2*rowsperpanel+1,] <- 10
for(i in 1:2){
for(j in 1:4){
lay.mat[(1+(i-1)*rowsperpanel):(i*rowsperpanel),(2+(j-1)*colsperpanel):(1+j*colsperpanel)] <- 4*(i-1) + j
}
}
layout(lay.mat)
growth.ests <- as.data.frame(xenic.strains)
colnames(growth.ests) <-'Strain'
growth.ests$a.linear <- NaN
growth.ests$a.linear.se <- NaN
growth.ests$a.all.linear <- NaN
growth.ests$a.all.linear.se <- NaN
growth.ests$a.x.linear <- NaN
growth.ests$a.x.linear.se <- NaN
growth.ests$a.r.linear <- NaN
growth.ests$a.r.linear.se <- NaN
growth.ests$a.sat <- NaN
growth.ests$a.sat.se <- NaN
growth.ests$h.sat <- NaN
growth.ests$h.sat.se <- NaN
for(i in 1:length(xenic.strains)){
subdat1 <- rbind(graz.data[graz.data$TargetOchConc!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
# Make the plot
plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.popn.24,na.rm=T)), ylim=c(min(subdat1$pred.mu.72.rel,na.rm=T),max(subdat1$pred.mu.72.rel,na.rm=T)),las=1)
abline(h=0,lty=1,col='gray50')
# Add labeling
if(i %in% c(1,2,5,6,7,8)){
text(x = 0-.06*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$pred.mu.72.rel,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)
}
if(i %in% c(3,4)){
text(x = 1.05*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$pred.mu.72.rel,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 2)
}
if(i == 8){
legend(x='topright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2)
}
if(i == 1){
mtext(expression(paste('Relativized ', italic('O. marina'),' growth rate (',d^-1,')')),side=2,cex=.8,at=-.2,line=2.5)
}
if(i == 6){
mtext(expression(paste('Prey availability (',italic('Ochromonas'),' · ',mL^-1, ')')),side=1,line=2.7,cex=.8,at=275000)
}
# Add xenic data
subdat.x <- subdat1[subdat1$Strain==xenic.strains[i] ,]
points(subdat.x$prey.popn.24,subdat.x$pred.mu.72.rel,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
# Add rice data
subdat.r <- subdat1[subdat1$Strain==rice.strains[i],]
points(subdat.r$prey.popn.24,subdat.r$pred.mu.72.rel, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
# Add linear model for subset of data
subdat1.trim <- subdat1[subdat1$prey.popn.24 < linear.fit.maxima[i],]
lm1 <- nls(pred.mu.72.rel ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
xset <- seq(from = 0, to = linear.fit.maxima[i])
yset <- summary(lm1)$parameters[1] * xset
#lines(xset,yset,lwd=1.5,col=col.strains[i])
growth.ests$a.linear[i] <- summary(lm1)$parameters[1]
growth.ests$a.linear.se[i] <- summary(lm1)$parameters[2]
# Add linear model for all data
subdat1.trim <- subdat1
lm1 <- nls(pred.mu.72.rel ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
xset <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T),length.out=100)
yset <- summary(lm1)$parameters[1] * xset
#lines(xset,yset,lwd=1.5,col='black',lty=2)
growth.ests$a.all.linear[i] <- summary(lm1)$parameters[1]
growth.ests$a.all.linear.se[i] <- summary(lm1)$parameters[2]
# Add linear model for xenic data
subdat1.trim <- subdat1[subdat1$Strain==xenic.strains[i],]
lm1 <- nls(pred.mu.72.rel ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
xset <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T),length.out=100)
yset <- summary(lm1)$parameters[1] * xset
lines(xset,yset,lwd=2.5,col=col.strains[i],lty=2)
growth.ests$a.x.linear[i] <- summary(lm1)$parameters[1]
growth.ests$a.x.linear.se[i] <- summary(lm1)$parameters[2]
# Add linear model for rice data
subdat1.trim <- subdat1[subdat1$Strain==rice.strains[i],]
lm1 <- nls(pred.mu.72.rel ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
xset <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T))
yset <- summary(lm1)$parameters[1] * xset
lines(xset,yset,lwd=2,col=col.strains[i],lty=1)
growth.ests$a.r.linear[i] <- summary(lm1)$parameters[1]
growth.ests$a.r.linear.se[i] <- summary(lm1)$parameters[2]
}
graz.data$Strain.numeric <- as.numeric(str_extract_all(graz.data$Strain, "\\d+"))
graz.data$Bact <- 'X'
graz.data[graz.data$Strain%in%c('584R','590R','1148R','1150R','1391R','1392R','1393R','2951R'),]$Bact <- "R"
graz.data$Strain.Bact <- paste(graz.data$Strain.numeric,graz.data$Bact,sep='.')
Barbaglia.data$Strain.Bact <- paste(Barbaglia.data$Strain,Barbaglia.data$Bact.short,sep='.')
# Append information for Ochromonas carbon content in pg C per cell
graz.data$C.perOch <- Barbaglia.data$C.perOch[match(graz.data$Strain.Bact,Barbaglia.data$Strain.Bact)]
# Convert ingestion rate from cells to C, divide by 1000 to get ng per day
graz.data$prey.i.24.C <- graz.data$prey.i.24*graz.data$C.perOch/1000
# Append information for Ochromonas carbon content in pg N per cell
graz.data$N.perOch <- Barbaglia.data$N.perOch[match(graz.data$Strain.Bact,Barbaglia.data$Strain.Bact)]
# Convert ingestion rate from cells to C, divide by 1000 to get ng per day
graz.data$prey.i.24.N <- graz.data$prey.i.24*graz.data$N.perOch/1000
par(mar=c(2,2.5,.7,.5))
rowsperpanel <- 5
colsperpanel <- 5
lay.mat <- matrix(data=NA,nrow=2*rowsperpanel+1,ncol=4*colsperpanel+1)
lay.mat[,1] <- 9
lay.mat[2*rowsperpanel+1,] <- 10
for(i in 1:2){
for(j in 1:4){
lay.mat[(1+(i-1)*rowsperpanel):(i*rowsperpanel),(2+(j-1)*colsperpanel):(1+j*colsperpanel)] <- 4*(i-1) + j
}
}
layout(lay.mat)
for(i in 1:length(xenic.strains)){
subdat1 <- rbind(graz.data[graz.data$TargetOchConc!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
# Make the plot
plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.i.24.C,na.rm=T)), ylim=c(min(subdat1$pred.mu.72.rel,na.rm=T),max(subdat1$pred.mu.72.rel,na.rm=T)),las=1)
abline(h=0,lty=1,col='gray50')
# Add xenic data
subdat.x <- subdat1[subdat1$Strain==xenic.strains[i] ,]
points(subdat.x$prey.i.24.C,subdat.x$pred.mu.72.rel,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
# Add rice data
subdat.r <- subdat1[subdat1$Strain==rice.strains[i],]
points(subdat.r$prey.i.24.C,subdat.r$pred.mu.72.rel, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
# Add linear model for xenic data
subdat1.trim <- subdat1[subdat1$Strain==xenic.strains[i],]
lm1 <- nls(pred.mu.72.rel ~ a*prey.i.24.C, data=subdat1.trim, start = list(a = 0.001))
xset <- seq(from = 0, to = max(subdat1.trim$prey.i.24.C,na.rm=T),length.out=100)
yset <- summary(lm1)$parameters[1] * xset
lines(xset,yset,lwd=2.5,col=col.strains[i],lty=2)
# Add linear model for rice data
subdat1.trim <- subdat1[subdat1$Strain==rice.strains[i],]
lm1 <- nls(pred.mu.72.rel ~ a*prey.i.24.C, data=subdat1.trim, start = list(a = 0.001))
xset <- seq(from = 0, to = max(subdat1.trim$prey.i.24.C,na.rm=T),length.out=100)
yset <- summary(lm1)$parameters[1] * xset
lines(xset,yset,lwd=2,col=col.strains[i],lty=1)
# Add labeling
if(i %in% c(1,5:8)){
text(x = 0-.06*max(subdat1$prey.i.24.C,na.rm=T), y = max(subdat1$pred.mu.72.rel,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)}
if(i %in% c(2:4)){
text(x = 1.06*max(subdat1$prey.i.24.C,na.rm=T), y = max(subdat1$pred.mu.72.rel,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 2)
}
if(i == 8){
legend(x='bottomright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2, lty = c(2,1))
}
if(i == 1){
mtext(expression(paste('Relativized ', italic('O. marina'),' growth rate (',d^-1,')')),side=2,cex=.8,at=-.18,line=2.5)
}
if(i == 6){
mtext(expression(paste('Ingestion rate (ng ',italic('Ochromonas'),' C · ',italic('O. marina')^-1,' · ',d^-1,')')),side=1,line=2.7,cex=.8,at=4.5)
}
}
par(mar=c(2,2.5,.7,.5))
rowsperpanel <- 5
colsperpanel <- 5
lay.mat <- matrix(data=NA,nrow=2*rowsperpanel+1,ncol=4*colsperpanel+1)
lay.mat[,1] <- 9
lay.mat[2*rowsperpanel+1,] <- 10
for(i in 1:2){
for(j in 1:4){
lay.mat[(1+(i-1)*rowsperpanel):(i*rowsperpanel),(2+(j-1)*colsperpanel):(1+j*colsperpanel)] <- 4*(i-1) + j
}
}
layout(lay.mat)
for(i in 1:length(xenic.strains)){
subdat1 <- rbind(graz.data[graz.data$TargetOchConc!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
# Make the plot
plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.i.24.N,na.rm=T)), ylim=c(min(subdat1$pred.mu.72.rel,na.rm=T),max(subdat1$pred.mu.72.rel,na.rm=T)),las=1)
abline(h=0,lty=1,col='black')
# Add xenic data
subdat.x <- subdat1[subdat1$Strain==xenic.strains[i] ,]
points(subdat.x$prey.i.24.N,subdat.x$pred.mu.72.rel,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
# Add rice data
subdat.r <- subdat1[subdat1$Strain==rice.strains[i],]
points(subdat.r$prey.i.24.N,subdat.r$pred.mu.72.rel, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
# Add linear model for xenic data
subdat1.trim <- subdat1[subdat1$Strain==xenic.strains[i],]
lm1 <- nls(pred.mu.72.rel ~ a*prey.i.24.N, data=subdat1.trim, start = list(a = 0.001))
xset <- seq(from = 0, to = max(subdat1.trim$prey.i.24.N,na.rm=T),length.out=100)
yset <- summary(lm1)$parameters[1] * xset
lines(xset,yset,lwd=2.5,col=col.strains[i],lty=2)
# Add linear model for rice data
subdat1.trim <- subdat1[subdat1$Strain==rice.strains[i],]
lm1 <- nls(pred.mu.72.rel ~ a*prey.i.24.N, data=subdat1.trim, start = list(a = 0.001))
xset <- seq(from = 0, to = max(subdat1.trim$prey.i.24.N,na.rm=T),length.out=100)
yset <- summary(lm1)$parameters[1] * xset
lines(xset,yset,lwd=2,col=col.strains[i],lty=1)
# Add labeling
if(i %in% c(1,5,6,8)){
text(x = 0-.06*max(subdat1$prey.i.24.N,na.rm=T), y = max(subdat1$pred.mu.72.rel,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)}
if(i %in% c(2:4,7)){
text(x = 1.06*max(subdat1$prey.i.24.N,na.rm=T), y = max(subdat1$pred.mu.72.rel,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 2)
}
if(i == 8){
legend(x='bottomright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2, lty = c(2,1))
}
if(i == 1){
mtext(expression(paste('Relativized ', italic('O. marina'),' growth rate (',d^-1,')')),side=2,cex=.8,at=-.18,line=2.5)
}
if(i == 6){
mtext(expression(paste('Ingestion rate (ng ',italic('Ochromonas'),' N · ',italic('O. marina')^-1,' · ',d^-1,')')),side=1,line=2.7,cex=.8,at=.9)
}
}
Barbaglia.data$Bact.short <- factor(Barbaglia.data$Bact.short,levels=c('X','R'))
# Useful features:
## C.from.grazing.perCperday = carbon obtained from grazing per Ochromonas carbon per day
## P_I.perC = carbon obtained from photosynthesis per Ochromonas carbon per day
Grazing
#Barbaglia.data$Strain.Bact.Light
par(mar=c(4,4,1,1))
attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'Best.a', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
plot(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,las=1, xlab='Ochromonas attack rate', ylab='O. marina attack rate',type='n',ylim=c(0.3,1.6),xlim=c(0.3e-6,4.7e-6))
arrows(attack.summ$Best.a-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$Best.a+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
subdat <- attack.summ[attack.summ$Strain==strains[i],]
lines(subdat$Best.a,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=c(24,21)[as.factor(attack.summ$Bact.short)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)
pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' clearance rate')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.95),xlim=c(0.3e-6,4.7e-6))
arrows(attack.summ$Best.a-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$Best.a+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
subdat <- attack.summ[attack.summ$Strain==strains[i],]
lines(subdat$Best.a,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)
legend(x='topright',inset=c(0.01,0.01),legend=xenic.strains,pch=pch.strains, cex=.8, pt.cex=cex.strains, col=col.strains, pt.bg=col.strains, title='Strain')
Chlorophyll
par(mar=c(4,4,1,1))
attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'chl.percell', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' chlorophyll per cell')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.8),xlim=c(0.12,0.34))
arrows(attack.summ$chl.percell-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$chl.percell+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
subdat <- attack.summ[attack.summ$Strain==strains[i],]
lines(subdat$chl.percell,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)
legend(x='topright',inset=c(0.01,0.01),legend=xenic.strains,pch=pch.strains, cex=.8, pt.cex=cex.strains, col=col.strains, pt.bg=col.strains, title='Strain')
Both together
par(mar=c(4,4,1,0.5),mfrow=c(1,2))
attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'Best.a', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' clearance rate')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.95),xlim=c(0.3e-6,4.7e-6))
arrows(attack.summ$Best.a-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$Best.a+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
subdat <- attack.summ[attack.summ$Strain==strains[i],]
lines(subdat$Best.a,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)
legend(x='topright',inset=c(0.01,0.01),legend=xenic.strains,pch=pch.strains, cex=.8, pt.cex=cex.strains, col=col.strains, pt.bg=col.strains, title='Strain')
par(mar=c(4,2.5,1,2))
attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'chl.percell', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' chlorophyll per cell')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.95)) #,xlim=c(0.12,0.34)
arrows(attack.summ$chl.percell-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$chl.percell+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
subdat <- attack.summ[attack.summ$Strain==strains[i],]
lines(subdat$chl.percell,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)
#legend(x='topright',inset=c(0.01,0.01),legend=xenic.strains,pch=pch.strains, cex=.8, pt.cex=cex.strains, col=col.strains, pt.bg=col.strains, title='Strain')
Total carbon yield
par(mar=c(4,4,1,1))
## C.from.grazing.perCperday = carbon obtained from grazing per Ochromonas carbon per day
## P_I.perC = carbon obtained from photosynthesis per Ochromonas carbon per day
Barbaglia.data$TotC.perC <- Barbaglia.data$C.from.grazing.perCperday + Barbaglia.data$P_I.perC
Barbaglia.data$TotC.percell <- Barbaglia.data$TotC.perC * Barbaglia.data$C.perOch
attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'TotC.percell', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' carbon yield rate')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.95),xlim=c(0,250))
arrows(attack.summ$TotC.percell-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$TotC.percell+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
subdat <- attack.summ[attack.summ$Strain==strains[i],]
lines(subdat$TotC.percell,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)
legend(x='topright',inset=c(0.01,0.01),legend=xenic.strains,pch=pch.strains, cex=.8, pt.cex=cex.strains, col=col.strains, pt.bg=col.strains, title='Strain')
Combining this with the attack rate figure
par(mar=c(4,4,1,0.5),mfrow=c(1,2))
attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'Best.a', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' clearance rate')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.95),xlim=c(0.3e-6,4.7e-6))
arrows(attack.summ$Best.a-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$Best.a+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
subdat <- attack.summ[attack.summ$Strain==strains[i],]
lines(subdat$Best.a,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)
legend(x='topright',inset=c(0.01,0.01),legend=xenic.strains,pch=pch.strains, cex=.8, pt.cex=cex.strains, col=col.strains, pt.bg=col.strains, title='Strain')
par(mar=c(4,2.5,1,2))
attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'TotC.percell', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' carbon yield rate')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.95),xlim=c(0,240))
arrows(attack.summ$TotC.percell-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$TotC.percell+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
subdat <- attack.summ[attack.summ$Strain==strains[i],]
lines(subdat$TotC.percell,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)
T-tests
# Make a holding table
t.test.results <- as.data.frame(xenic.strains)
colnames(t.test.results) <- "Strain"
t.test.results$Och.graz.tval <- NaN
t.test.results$Och.graz.pval <- NaN
t.test.results$Oxy.graz.numSEs <- NaN
for(i in 1:length(xenic.strains)){ # Working strain by strain
strainchoice <- xenic.strains[i]
# Gina's data: Is there a difference in clearance rate between xenic and rice data?
# Negative t-value for Gina's data means the rice strains had *higher* clearance rates than the xenic ones
gina.subdat <- Barbaglia.data[Barbaglia.data$Light==100 & Barbaglia.data$Strain==strainchoice,]
t.test.results$Och.graz.tval[i] <- t.test(Best.a ~ Bact.short,data=gina.subdat)$statistic
t.test.results$Och.graz.pval[i] <- t.test(Best.a ~ Bact.short,data=gina.subdat)$p.value
}
# AG's data: Is there a difference in clearance rate between xenic and rice data?
# Negative test value means rice strains had *higher* clearance rates
subdat <- as.data.frame(cbind(attack.ests$Strain,attack.ests$a.r.linear,attack.ests$a.r.linear.se))
colnames(subdat) <- c('Strain','attack','se')
subdat$Bact <- 'R'
subdat2 <- as.data.frame(cbind(attack.ests$Strain,attack.ests$a.x.linear,attack.ests$a.x.linear.se))
colnames(subdat2) <- c('Strain','attack','se')
subdat2$Bact <- 'X'
subdat <- rbind(subdat2,subdat)
subdat <- subdat[order(subdat$Strain),]
for(i in 1:length(xenic.strains)){
subdat2 <- subdat[subdat$Strain==xenic.strains[i],]
max.se <- max(subdat2$se)
delta.attack <- subdat2[subdat2$Bact=='X',]$attack - subdat2[subdat2$Bact=='R',]$attack
t.test.results$Oxy.graz.numSEs[i] <- delta.attack / max.se
}
print(t.test.results)
## Strain Och.graz.tval Och.graz.pval Oxy.graz.numSEs
## 1 584 3.0993720 0.0647563383 -8.778937
## 2 590 -0.2824363 0.7981552021 -2.951987
## 3 1148 -0.1434081 0.8937054002 1.835591
## 4 1150 0.6158596 0.5903486015 1.268617
## 5 1391 11.5989292 0.0003357685 -10.426857
## 6 1392 5.4398070 0.0075964682 -5.033688
## 7 1393 9.8461959 0.0006750262 -3.972854
## 8 2951 18.9398352 0.0009854816 -15.845488